Nonlinear elasticity of monolayer graphene 
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By combining continuum elasticity theory and tight-binding atomistic simulations, we work out 
the constitutive nonlinear stress-strain relation for graphene stretching elasticity and we calculate 
all the corresponding nonlinear elastic moduli. Present results represent a robust picture on elastic 
behavior of one-atom thick carbon sheets and provide the proper interpretation of recent experi- 
ments. In particular, we discuss the physical meaning of the effective nonlinear elastic modulus there 
introduced and we predict its value in good agreement with available data. Finally, a hyperelastic 
softening behavior is observed and discussed, so determining the failure properties of graphene. 

PACS numbers: 62.25.-g, 62.20.D-, 46.70.Hg 



The elastic properties of graphene have been recently 
determined by atomic force microscope nanoindentation 
[U |2], measuring the deformation of a free-standing 
monolayer as sketched in Figjl] (top). In particular, in 
Ref.[l] the experimental force-deformation relation has 
been expressed as a phenomenological nonlinear scalar 
relation between the applied stress (a) and the observed 
strain (e) 
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where E and D are, respectively, the Young modulus 
and an effective nonlinear (third-order) elastic modu- 
lus of the two dimensional carbon sheet. The reported 
experimental values are: E = 340 ± 40 Nm~^ and 
D = — 690±120Nm~^. While the first result is consistent 
with previous existing data UHl], the above value for D 
represents so far the only available information about the 
nonlinear elasticity of a one- atom thick carbon sheet. 

Although nonlinear features are summarized in Eq.Q 
by one effective parameter continuum elasticity theory 
predicts the existence of three independent third-order 
parameters Cijk for graphene, as reported below. In 
other words, while Eq.Q represents a valuable effective 
relation for the interpretation of a complex experiment 
[1], it must be worked out a more rigorous theoretical 
picture in order to properly define all the nonlinear elas- 
tic constants of graphene and to understand the physical 
meaning of D. This corresponds to the content of the 
present Letter where we investigate the constitutive non- 
linear stress-strain relation of graphene, by combining 
continuum elasticity and tight-binding atomistic simula- 
tion (TB-AS) |8^. 

To obtain the nonlinear stress-strain relation of an elas- 
tic membrane, we need at first to elaborate an expression 
for the corresponding strain energy function U (per unit 
area). Since, as illustrated in Fig. 1 (bottom), the under- 
lying lattice is hexagonal, it is useful to consider the co- 
ordinate set a = X -\- iy and ^ = x — iy [9^^ where the 
X and y directions are respectively identified with the 



zig-zag (zz) and the armchair (ac) directions. Because 
the strain energy function is invariant under a rotation 
of 7r/3 about the z-axis (normal to the suspended mono- 
layer), there are two linear moduli (the two-dimensional 
Young modulus E and Poisson ratio u) and three non- 
linear independent elastic coefficients (A^, i = 1,2,3) all 
expressed in units of force/length; we easily proved that 
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Figure 1: (color online) Top: schematic representation of the 
indentation of a suspended monolayer graphene (side view). 
Bottom: definition of zig-zag and armchair directions (top 
view) . The gray-scale shading of the monolayer graphene pic- 
torially represents the radially symmetric strain field gener- 
ated by indentation. 
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where 6^/3 = 6^ 
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2iexy, and = 



^xx — ^yy — 2^63;^. In order to further proceed we must 
better focus the strain definition which in elasticity the- 
ory is twofold: we can introduce the so-called small strain 
tensor e = ^{Vu-\- Vil^), being u the displacement field, 
or the Lagrangian strain i) = ^(Vu -\- Vu^ + ViTVu). 
While e takes into account only the physical nonlinearity 
features (it describes a nonlinear stress-strain dependence 
observed in regime of small deformation), fj describes any 
possible source of nonlinearity, i.e. it includes both phys- 
ical and geometrical (large deformation) ones. 

We start using e in Eq.([2| and we get the nonlinear 
elastic coefficients Ai, A2 and A3 which are related to 
the third order elastic constants Cm, C222 and C112, as 
customarily defined in crystal elasticity (TOj , through the 
following relations 



C222), A2 — -(C222 — C112), 



A3 — — (2C111 — C222 + 3C112). 

The strain energy function is finally obtained as 
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where we set Caa^fBfB = Tr (e^) and e^^ = (Tre)^. The 
stress-strain nonlinear constitutive equation for in-plane 
stretching is straightforwardly obtained by T = dU/de^ 
where T is the Cauchy stress tensor. 

Since the analysis of the experimental data provided 
m Ref.[l] through Eq.Q is assuming an applied uniax- 
ial stress, we now suppose to apply a uniaxial tension 
afi along the arbitrary direction ft = cosOcx + sin6>e^, 
where Cx and Cy are the unit vectors along the zig-zag 
and the armchair directions, respectively (see Figjl], bot- 
tom). Under this assumption we get: T = af^n (g) n, 
with in-plane components defined as T^x = cos^ ^, 
Txy = cr^cos^sin^, and Tyy = afisin^O. Similarly, by 
inverting the nonlinear constitutive equation we find the 
corresponding strain tensor and the relative variation of 
length Cfi = n ■ en along the direction n. By combining 
these results, we obtain the stress-strain relation along 
the arbitrary direction n 



Duel 
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where Dfi is given by 



Dn = l{l- A3 + ^ (1 - ^) (1 + vf A2 (6) 
+3 (2 cos^ (9-1) (16 cos"^ 19-16 cos^ i9 + l) (1 + uf Ai 



If we set fi = Cx (i.e. = 0), we get the nonlinear 
modulus D^^^) for stretching along the zig-zag direction 



D 
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Similarly, by setting ft = Cy (i.e. ^ = 7r/2), we obtain the 
nonlinear modulus for stretching along the arm- 

chair direction 

D('^'=)=De-, = -3(l + /.f Ai + 5(1-j.)(1 + j.)2a2 



A3 
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We observe that the above expression for D'^^^^ apply for 
all stretching directions defined by the angles 6 = kir/S 
{keZ), while D^^^^ holds for the angles (9 = 7r/6 + /c7r/3. 

Since the nanoindentation experiments generate a 
strain field with radial symmetry [T, as sketched in 
Figjljbottom), in order to get the unique 5ca/ar nonlinear 
elastic modulus appearing in Eq.(l) we need to average 
the expression of Dfi over 0. This procedure leads to 

1 r^TT nizz) I r)(ac) 
(D,) = ^J^ D,,0= ± 

= 5 (1 - [(l+,/)»A2 + Aj] (9) 

proving that the experimentally determined nonlinear 
modulus actually corresponds to the average value of the 
moduli for the zig-zag and armchair directions. 

We now repeat the above procedure by using the La- 
grangian strain 7): even in this case it is demonstrated 
that the strain energy function is given by the very same 
Eq.Q where e is replaced by fj and the Cijk by the La- 
grangian third-order moduli C^-^. By imposing the iden- 
tity /7(e) = Uifj) (where the Lagrangian strain can be 
written in term of the small strain by 7) = e + ^e^ [EH US]) 



we obtain the conversion rules: = Cm — jz^? 

^2^22 = ^^222 — 13^2 9 C112 = ^^112 — = Dfi — 

(for any n) and (D^) = (Dfi) - \E. The constitu- 
tive equation can be finally derived in the form = 
dU/dfj^ where is the second Piola-Kirchhoff stress 
tensor. Hereafter we will refer to the small strain and 
Lagrangian scalar nonlinear modulus by {Dfi) and {D^)^ 
respectively. They both will be compared with the ex- 
perimental parameter D of Eq.Q. The analysis below 
will identify the actual theoretical counterpart of D. 

The important result summarized in Eq.(|9| (as well as 
in its Lagrangian version) implies that the scalar non- 
linear modulus can be obtained by the third-order elas- 
tic constants (as well as the linear ones). They can 
be computed through the energy-vs-strain curves corre- 
sponding to suitable homogeneous in-plane deformations, 
thus avoiding a technically complicated simulation of the 
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Table I: Relationship am ong the energy expansion coefficients 
[7*^^-* and U^^^ of Eq.(lO) and the elastic moduli of graphene 



for four in-plane deformations (see text). 
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Table II: Small strain and Lagrangian nonlinear elastic mod- 
uli of graphene in units of Nm~^. 



Small strain 



Lagrangian 



Cm 


-1689.2 


Cm 


-2724.7 


C222 


-1487.7 


^222 


-2523.2 


C112 


-484.1 


^112 


-591.1 




-696.2 


jjCizz) 


-1163.7 




-469.6 




-937.9 



nanoindentation experiment. Therefore, the following in- 
plane deformation have been applied: (i) an uniaxial de- 
formation ( along the zig-zag direction, corresponding 
to a strain tensor e-^^^ = (dixSjx; (ii) an uniaxial de- 
formation C along the armchair direction, corresponding 
to a strain tensor e^^^^ = (SiySjy; (iii) an hydrostatic 
planar deformation corresponding to the strain tensor 
e-J^ = (5ij; (iv) a shear deformation corresponding to 

an in-plain strain tensor e-^^ = ( {SixSjy + SiySjx)- 

In this work the needed energy-vs-strain curves have 
been determined by TB-AS, making use of the tight- 
binding representation by Xu et al. [13 . A periodically 
repeated square cell containing 400 carbon atoms was de- 
formed as above. For any given applied deformation, full 
relaxation of the internal degrees of freedom of the sim- 
ulation cell was performed by zero temperature damped 
dynamics until interatomic forces resulted not larger than 
0.5 • 10~^^eV/A. Similar calculations were repeated by 
using a smaller simulation cell containing 200 atoms and 
by relaxing the system through simulated annealing. No 
deviation from data here reported were observed. 

For the deformations e,- , e,- , e,-^^ and e,-^^ the elastic 
energy of strained graphene can be written in terms of 
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Figure 2: (color online) Strain energy density U, obtained by 
TB-AS, as function of the strain parameter ( corresponding 
to the four homogeneous deformations summarized in Tab|l] 
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where Uq is the energy of the unstrained configuration. 
Since the expansion coefficients [/^^^ and U'^^^ are related 
to elastic moduli as summarized in Tabjlj a straightfor- 
ward fit of Eq.( 10 ) has provided the full set of linear mod- 



uli and third order elastic constants, while the shear de- 
formation was used to confirm the isotropy of the lattice 
in the linear approximation. Each energy-vs-strain curve, 
shown in Fig|2j has been computed by TB-AS as above 
described, by increasing the magnitude of ( in steps of 
0.005 up to a maximum strain |Cmax| = 0.055. Arrows 
in Figj2] indicate the different nonlinear behavior along 
the zz and ac directions. A similar fitting procedure was 
carried out by computing all the components of the stress 
tensor for the homogeneous deformations, obtaining no 
quantitative difference in the calculated moduli. 

The outputs of the fitting procedure are reported in 
Tab nil where the full set of third order elastic con- 
stants of monolayer graphene is shown. We remark that 
Cm is different than C2225 i-e a monolayer graphene is 
isotropic in the linear elasticity approximation, while it 
is anisotropic when nonlinear features are taken into ac- 
count. By inserting the elastic constants Cijk of Tab|TI| 
into Eqs.(|3|, ^ and ([8|, we also obtained the nonlinear 
moduli for both the zz and ac directions. 

In Tab|III| we report the values of the calculated elas- 



Table III: Linear and nonlinear elastic moduli of graphene in 
units of Nm~^ (we remark that the Poisson ratio v is dimen- 
sionless) . 
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Figure 3: (Color online) Theoretical (present work) and 
experimental (see Ref.^j) stress-strain curves, as defined in 
Eq.g. Shaded area represents the experimental error. 



in the present problem. In addition, by means of Fig j3]we 
can determine the failure stress (maximum of the stress- 
strain curve) (7/ = —E^ I {^l\D^\ corresponding to a pre- 
dicted failure stress as high as 42.4 Nm~^. This result is 
in excellent agreement with the experimental value 42 ±4 
Nm~^, reported in Ref.[l]. These values correspond to 
the failure strength of a two-dimensional system. In order 
to draw a comparison with bulk materials, we can define 
an effective three-dimensional failure stress cr^^ = crj/d, 
where d can be taken as the interlayer spacing in graphite. 
By considering d = 0.335 nm [20 , we obtain <jj^ = 130 
GPa. This very high value, exceeding that of most ma- 
terials (even including other stiff carbon-based systems, 
e.g. multi- walled nanotubes [211), motivates the use of 
one-atom thick carbon layers as possible reinforcement 
in advanced composites. 

We acknowledge financial support by the project 
MIUR-PON "CyberSar". 



tic moduli, together with the available experimental and 
theoretical data. The present TB-AS value for E is in 
reasonable agreement with literature ^ 13 [H [15] , while 
the value of v is larger than most of the ab-initio results 
la HI HSl HI] (but for the result in Ref. [16]). While 
this disagreement is clearly due to the empirical charac- 
ter of the adopted TB model (where, however, no elastic 
data were inserted in the fitting data base), we remark 
that the values of (I^fi) and {E)i) predicted by means of 
Eq.([9| are affected by only 10% by varying v among the 
values shown in Tab lllll 

Tab |III| shows that the predicted {Pri) is much closer to 
the experimental value D than its Lagrangian counter- 
part {D^). This seems to suggest that, measurements 
in Ref.[l were performed in the physical nonlinearity 
regime (small strain formalism), rather than in the ge- 
ometrical nonlinearity one (Lagrangian formalism), as 
also confirmed by the excellent agreement shown in Fig|3] 
commented below. We further observe that the nega- 
tive sign of all the nonlinear elastic moduli proves that 
graphene is an hyperelastic softening system (i.e. D <{)). 
Therefore, as recently established [18, 19^, the present 
nonlinear model plays a crucial role in determining the 
failure behavior of the graphene membrane. 

In order to substantiate the above statement, we show 
in FigjSj the graphene stress-strain curve, as defined in 
Eq.Q. Both the theoretical and experimental curves 
have been obtained by using the Young modulus and the 



scalar nonlinear coefficient as reported in Tab III We re- 
mark that in Fig js] the small strain {Dfi) value was used. 
The agreement between the experimental curve and the 
theoretical (small strain) one is remarkable. This con- 
firms that likely only physical nonlinearities are at work 
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